How Multivalency controls Ionic Criticality 
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To understand how multivalency influences the reduced critical temperatures, T* (z), and densities, 
p*c{z), of 2 : 1 ionic fluids, we study equisized hard-sphere models with z—l — i. Following Debye, 
Hiickel and Bjerrum, association into ion clusters is treated with, also, ionic solvation and excluded 
volume. In good accord with simulations but contradicting integral-equation and fleld theories, T* 
falls when z increases while p* rises steeply: that 80 — 90% of the ions are bound in clusters near 
Tc serves to explain these trends. For z^l interphase Galvani potentials arise and are evaluated. 
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Multivalent ions play a significant role in condensed- 
matter, physicochemical, biophysical and, via the plasma 
transition, astrophysical contexts The effects of mul- 
tivalency are, however, often hard to comprehend. One 
central issue — relevant to electrolyte solutions, molten 
salts, Hquid metals, and dense plasmas Q — arises in 
Coulomb-driven phase separation. The most basic model 
for such ionic fluids consists of TV = pV hard-core spher- 
ical ions of various species cr in a volume V of uniform 
dielectric constant D, with — Pa-V ions of diameter 
flcr carrying charges qa = Zaqa, where go is an elementary 
charge. In the simple equisized z:l charge asymmetric 
primitive models (C^APMs), on which we focus here, one 
has a = -\-,—, a+ = a_, and q+ — zqo, q- — —qo- The ba- 
sic energy scale and associated reduced temperature and 
density are then e — zq^/ Da, T* = kBT / e, p* =pa^. 

Monte Carlo simulations |2l] show that (at least for 
z<5) the C^APMs exhibit "gas-Hquid" phase separation; 
furthermore, the critical parameters, T*{z) and Pc{z), 
are found to reasonable precision : see Table and the 
open circles in Figs. ^ and El One observes that T*{z) 
falls with increasing z, while p*{z) rises sharply. But 
we ask : How can these trends be understood? Or ac- 
counted for semiquantitatively? To address this issue we 
review briefly previous work, including a pioneering fleld- 
theoretic attack 0, and then report on a recent study 
01 which we believe provides signiflcant insight. This ex- 
tends an earlier analysis I 0| for the symmetric z — 1 re- 
stricted primitive model (RPM) that was founded on the 
original Debye-Hiickel (DH) approach but incorporated 
(i) Bjerrum ion pairs and (ii) their solvation in the resid- 
ual ionic fluid. For z = 2 and 3 larger ion clusters, trimers 
and tetramers, must be included 0, 0| ; but then explicit 
results are also obtained for the interphase Galvani po- 
tential f3| that appears in any two-phase nonsymmetric 
ionic system [j,L3|. 

The field-theoretic analysis of Netz and Orland (NO) 
Q was designed to address z:l ionic fluids and colloids 
(z ^ 1) and to include correlations in a systematic man- 
ner. The Coulomb interaction, q^qr/r, was transformed 
to yield a functional integral over an auxiliary poten- 
tial (/'(r). At the {(fP') level the DH effective interaction. 
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FIG. 1: Reduced critical temperatures for z:l charge 
asymmetric equisized hard-core primitive model electrolytes 
(CzAPMs ) according to Monte Carlo (MC) simulations 2] ; 
Debye-Hiickel (DH) theory; fleld-theoretic approaches : NO 
i3j (with a factor j^) and NMF approximate integral equa- 
tions : MSA Q, SPB, and MPS |10l|: and the present DHBjCI 
and DHBjCIHC solvated ion-cluster theory |j|; See text. 



TABLE I: Monte Carlo (MC) estimates I3| for the reduced 
critical parameters for z:l equisized hard-sphere electrolytes, 
values calculated from DHBjCIHC theory (CI) H, and ap- 
proximate estimates based on ion cluster statistics : see text. 



critical temp. 
z MC CI 

1 4.933 5.569 

2 4.70 4.9O7 

3 4.10 4.334 



io'rc(2) 

Edh Emc 

5.45 4.935 
5.11 4.65 
4.85 4.44 



critical density 

MC CI 
7.50 2.6I4 
9.3 6.26i 

12.5 11.90 



lOVc(^) 
Ep E« 
2.72 2.37 
4.27 3.49 
6.96 5.40 



i^oH oc e "''/r is captured with 

K\T-{p„})^Mql/DksT)Y, ^Ip^- (1) 

The reduced free energy density, /(T; p) = —FlVksT, 
was computed to eighth order in (/) but a momentum cut- 
off is essential: NO adopted |kA| = 27r/a thereby incorpo- 
rating the ionic diameter and, for the z:l case, leading to 




FIG. 2: Reduced critical densities pl(z), for the Cz APM elec- 
trolyte as in Fig. U (except that the NO plot is not rescaled). 

k2q2 _ 47rp*/T*. Since this treatment of the hard cores 
is approximate, accurate predictions for T*{z) and p*{z) 
are not expected. Nevertheless, one might anticipate reli- 
able trends when z varies in contrast to DH theory which 
yields no dependence on z with (after I) 

DH : Kaa = 1, T* ^ ^, = l/647r ~ 0.005 . (2) 

In fact, as NO report, "the [predicted] deviations from DH 
theory are pronounced" for z > 1 : see the bold dashed 
plots in Figs, n] and 12 

But evidently the NO results are not merely quantita- 
tively wrong; the trends are quite incorrect since T* is 
asserted to rise rapidly (instead of falling) while p* falls 
sharply for small z—1 (instead of rising) and then in- 
creases but much too slowly. While one may blame the 
approximate treatment of the hard cores, we believe this 
is not the primary culprit. Indeed, a recent field-theoretic 
analysis paid closer attention to the ion- ion repulsions 0|; 
but the subsequent "new mean-field" (NMF) results still 
exhibit strong increases in T* and an overly weak varia- 
tion of p*: see the NMF plots in Figs. [Hand El 8]. 

Integral equation theories are hardly better : see 
Figs.[l]and|2l The mean spherical approximation (MSA), 
like DH theory, predicts no variation of T* and p* with z 
0. A symmetric Poisson-Boltzmann (SPB) theory 
does predict the correct falling and rising trends for T* 
and p*, but the degree of variation is woefully inadequate. 
Moreover, the modified Poisson-Boltzmann (MPB) ap- 
proximation, that the same authors ,10] argue should be 
more reliable, yields the wrong trend for T*. 

In order to better understand the effects of multiva- 
lency we turn to recent calculations 0, based on the 
solvated ion-cluster view of the C^APM near criti- 
cality that is supported 'pictorially' by simulations 0- 
In brief, the aim is to construct the free energy density, 
/(T; {pa-}), for ionic species a consisting (i) of + and — 
monomers, i.e., isolated, n+ =n_ = 1 single, unassociated 



FIG. 3: Coexistence curves predicted for z:l equisized prim- 
itive models by the DHBjCI and DHBjCIHC theories (solid 
and dashed lines, respectively) together with Monte Carlo es- 
timates (MC) based on [Tl| . 



ions of valency z+ = z and Z- = —1 ; (ii) a set of associ- 
ated primary clusters, a = 2, 3, . . ., dimers, trimers, . . ., 
each consisting of one "central" -I- ion and = cr — 1 
"satellite" counterions for a total of ria- — + 1 ions in 
a cluster of valency z^ — z — ma-; up to (iii), the largest 
primary cluster, the neutral or 'molecular' {z + l)-mer of 
one z+ ion and z negative ions 

For each species, / contains an ideal-gas term 
P'^{T,pcr), and an electrostatic term /J^'(r, {p^}), that, 
following DH, incorporates Cluster solvation in the par- 
tially associated Ionic fluid: this description is thus 
dubbed "DHBjCI" 3|. By adding a Uard Core free- 
volume term, /""^({pcr}), as in I, one may also account for 
those excluded volume effects not already encompassed 
in the basic solvation and association calculations 0,0]; 
so generating a "DHBjCIHC" theory (The effective 
HC virial coefficient B]^"'' = has been adopted 

[3,IB|-) Examination of Figs. and [2 reveals that these 
solvated ion-cluster theories are surprisingly successful! 
Not only are both the downward trend in T*{z) and the 
rapid rise of p^iz) well captured, but the quantitative 
agreement with each of the MC estimates is significantly 
better than achieved by other approaches. 

One must recognize that (all) these theories are of 
mean-field character: thus 5 to 15% over-estimates of 
T*{z) are to be expected. Indeed, neglected fiuctuations 
typically depress Tc by such amounts and also fiatten the 
coexistence curves as seen in Fig. |2l Second, note that 
the hard-core terms have a small effect on p*{z) while re- 
ducing T*{z) values by only 5—10%. Nevertheless, Fig.EI 
reveals that the liquid phases, especially for p* > 0.15, 
are sensitive to /^^: but, recall the discussion in I. In 
fact, the crucial feature of DHBj-type theories — not rep- 
resented in field-theoretic or standard integral-equation 
treatments — is the chemical equilibrium maintained be- 
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tween the cluster species via the Law of Mass Action: 

P<T = Km AT) P+ exp [pl""^ + mpL''_^ - ^C] , (3) 

for a = m + 1 > 2, with the excess chemical potentials 
/i^"' — —{d/dpcr)[f"^' +Y^afa']y while the association 
constants are taken as 

in which Sm,z{{T^i}) is the electrostatic energy of an iso- 
lated (m + l)-mer with satellite coordinates {r^}. The 
lower limits a and the condition £m^z — +oo for jr^ — | < 
a, represent hard cores. Following Bjerrum 0, the nec- 
essary cut-off radius R is chosen so that {dKm,z/dR) is 
minimal. The resulting 3-fold K2,z integral is managable 
but the 6-fold integral for if 3, 3 requires a Pade approx- 
imant study of the low-T expansion cross-checked to a 
part in 10^ by MC evaluations It transpires, how- 
ever, that T* and p* are insensitive to the Km z values 

a. 

Lastly, one needs to account for the solvation of all the 
ion species, a, by the free ions and charged clusters via 
the electrostatic terms [4, ^ ^] 

/-(T;{M) = ^E^E(IQ-P), (5) 

/— rn— — l 

where the tPnix) are related to the spherical Bessel func- 
tions ki{x) yi; the second sum requires the cluster electric 
multipole moments, Qfm' thermally averaged 0| over the 
ionic configurations that already enter in the Km,z{T). 

Finally, Oo- is an effective cluster diameter, i.e. the ra- 
dius of the approximating sphere (centered to minimize 
f^') that substitutes for the true, thermally fluctuating, 
hard-core exclusion domain: see I and y. One con- 
cludes, as in I, that a most reasonable choice for is the 
average over solid angle of the radial distance to the true 
exclusion surface of the ground-state cluster : this yields 
02 = (I -f I ln3 ~ 1.162)o, as = 1.250a and 04 = 1.375a. 
For z = l the values of T* and p* vary by less than ±2% 
over plausible alternatives for 02 Q; but the sensitivity 
to 03 and 04 for z = 2 and 3 is greater. As a result, this 
hard-to-avoid approximation contributes significantly to 
the overall quantitative uncertainties. 

From the total free energy /(T, {po-}), all thermody- 
namic properties follow One may then conclude 
from Figs. Hand El that the principal defect of the field- 
theoretic and integral-equation approaches is a failure to 
account effectively for strong ionic association near criti- 
cality. But can the actual trends of T* and p* with z be 
demonstrated in a direct, transparent way? To answer, 
consider the fractions, jja ~ ria-Nc^/N, of ions bound in 
clusters of n^r ions with pa- — (ycr/^cr)p- The critical point 



TABLE IL Inverse screening length k and fractions, j/c — 
UaNa/N, of ions in clusters of n„ ions at criticality, as per- 
centages, according to DHBjCIHC theory [J. 
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KcO, 




y- 


y2 


2/3 




1 


L04 


9.14 


9.14 


81.72 






2 


1.37 


1.31 


10.33 


15.43 


72.93 




3 


1.57 


0.34 


8.04 


3.32 


11.13 


77.17 



values that result from DHBjCIHC theory are dis- 

played in Tableim A significant fact is the rapid decrease 
in y'^, the fraction of unassociated 2+ ions, from 9.1 to 
1.3 to 0.3%. But more can be learned! 

To understand the variation of T*(z) let us regard the 
electrolyte in the critical region as a mixture of clusters 
with fixed mole fractions Xcr = {yalna)/^^{yT/'nr). A 
pair (cr, r) will either mutually repel or attract with pair- 
wise binding energies, say. Ear- Thus unlike monomers 
attract with e±=£. However, a dimer attracts only neg- 
ative monomers with £2- — {z — ^)£/z; but repels all 
z+ > +2 ions. Two dimers repel when z > 3; but one has 
£2,2/ £ — 0.586 and 0.345 for z = l and 2. And so on. 

To estimate T* for this mixture we adopt a van-der- 
Waals approach as in [2(d)]. Thus, for the overall cluster 
density p (= pY.aV<y/'^<y)-> take p/pkeT ~ Z{B^p) + 
Bi(T*)p with Z{u) = l+it+. . . in which the second virial 
coefficient has been decomposed as B(T*) = Bq + Bi{T*) 
where Bq (= 6oa^, say) represents the hard-core repul- 
sions while Bi{T*) embodies the attractions. Solving 
dpp = dpP = 0, as usual, yields p* and B* = Bi{T*) /boa^ . 
At low T, which is relevant here, one has 

Bi{T*) « - V b^ra^'xaXreMCr/T*) , (6) 

where e*^ = ecrr/e, while b„ra^ specifies the volume of 
mutual attractions: this vanishes if a and t repel. 

Now, the x+X- term dominates in Bi{T*) at low T 
with corrections of relative order {x\/ x+x^^^'^'^'^^^'^ 
for z = 1 and 2(a;2/a;+)e~^/^^"^ for z > 2. We may 
then calibrate Bi{T*)/a^ by using pure DH theory ^ for 
which, since association is not considered, .t+ = .x_ = |. 
Thereby we obtain the Edh estimates 

T;(z) ~ 1/(16 + I \n4xl{z)x'L{z)\) , (7) 

in which x'^ oc y'^ and xt oc yl follow from Table |n| 

The resulting predictions are listed in Table ID under 
Edh- In light of the heuristic nature of the arguments, 
they refiect the trend of the MC and CI values surpriz- 
ingly well. Certainly the contention that association is a 
prime factor is well confirmed. By replacing 16 by 20.27 
(or 17.96) in Q, and the factor 4 by l/x+(l)a;i(l), one 
caHbrates Bi(T*) on the MC (or CI) values for the RPM. 
Column Emc in Table ID lists the MC-calibrated values : 
for z = 2 and 3 these match the Monte Carlo estimates 
to within 1% and 8%, respectively (l^ . 
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FIG. 4: Reduced Galvani potentials, A(t> = qolS.(t>/kBT ^ vs. 
T/Tc for z:l electrolytes according to pure DH theory (dotted) 
and DHBjCI(HC) theories : solid (dashed) plots. 



Specifically, we have summarized briefly analytical cal- 
culations for 3:1, 2:1 and 1:1 equisized charged hard- 
sphere fluids [1| that, for the flrst time, reasonably reflect 
the true variation of critical temperatures and densities, 
T*{z) and p*c(z) (as revealed by simulations Q). On 
that basis, supported by analysis that correlates T*{z) 
and p*c^{z) with the increasingly depleted populations of 
free ions and charged clusters as z increases, it is clear 
that recognizing ionic association is inescapable for a suc- 
cessful theory. Previous treatments |j,la, lEUil, lacking 
allowance for ion clusters fail seriously. The ion-cluster 
solvation theories also yield quantitative results for the 
interphase Galvani potentials. 

National Science Foundation support via Grants CHE 
99-81772 and 03-01101 is gratefully acknowledged. 



Now, for the critical density, the signiflcance of ion 
pairing is already clear in pure DHBj theory for the RPM 
p. The heavy depletion of the free ions (which, in DHBj 
theory, drive the transition alone) means that to reach 
criticality the overall density p{=p+-^p- + 2p2) must be 
increased until the DH criterion p*^ + p*_= p'^% = 1 /GAtt 
is met : see Does the same depletion-by-association 
mechanism account for the z-dependence of Pci^)? 

To progress, rewrite 1^ generally as k'^o^ — Anp^ /T* , 
with the effective, depleted ionic density 

P^=P*Y] zlya(T;{p„})/zna . (8) 

If one accepts the DH criterion and uses Table ^ the 
estimates Ep, in TableHl result. Although these fall short 
of the Monte Carlo values by 74, 54 and 44% for z = 
1—3, they reproduce the accelerating increase with z (by 
factors 1.57, 1.63 vs. 1.24, 1.34). 

An alternative approach adopts the DH value Kca=l: 
see 10 but note from Table |ll| that DHBjCIHC theory 
implies that rises from 1.04 for the RPM to 1.57 
for z = 3. Then using the Edh values for T* , in Table 
leads to the E^ predictions for pl{z) : these are all rather 
low but the increases with z, by factors 1.47, 1.55, again 
reflect the correct behavior. 

Finally, we note that the Galvani potential, A0, that 
arises between coexisting phases in charge asymmetric 
fluids is readily calculated 0, Q| . The predictions from 
pure DH theory are shown dotted in Fig. : one flnds 
A0DH oc (1 — z~^). The other plots result from the DH- 
Bj CI and DHBjCIHC theories Surprizingly, the cal- 
culations suggest no clear trend with z. It is natural to 
conjecture that A(/) vanishes as Go{Tc — moreover 
to the extent that the expected mean-fleld value f3—\\s 
realized, the present results support this. 

In conclusion, we have elucidated the mechanisms un- 
derlying how multivalency influences critical behavior. 
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